#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"

extern handle_triple();

logical rayleigh_taylor()
{
return ((fluid_select == 24) || (fluid_select == 25));
}

logical rise_bubble()
{
return ((fluid_select==21) || (fluid_select==41) || (fluid_select==42));
}

logical use_square_box()
{
return ((fluid_select==21)&&(t_stop<10.0));
}

logical bubble_merge()
{
return (fluid_select == 27);
}

logical drop_merge()
{
return ((fluid_select == 60) || (fluid_select==160) || (fluid_select==161));
}

logical drop_periodic()
{
return (fluid_select == 61);
}

logical drop_jet()
{
return (fluid_select == 62);
}

logical oil_merge()
{
return ((fluid_select == 28) || (fluid_select == 29) || 
        (fluid_select==40) || (fluid_select == 50));
}

logical stratified_flow()
{
return (fluid_select==170);
}

logical dambreak()
{
return (fluid_select==90);
}

logical periodic_test()
{
return (fluid_select == 70);
}

logical three_phase()
{
return ((fluid_select==160)||(fluid_select==161)||(oil_merge())||
        (fluid_select==147) || (stratified_flow()) );
}

logical poolflag()
{
logical temp;

  temp=FALSE;
  if ((mx_type==REFLECT_BDRY)&&(my_type==REFLECT_BDRY))
   {
    if ((fluid_select==43)||(fluid_select==47)||(dambreak()) )
     temp=TRUE;
   }
  return temp;
}

logical feed_back()
{
return (fluid_select != 161);
}

logical use_air_water()
{
return ((fluid_select == 21) || (fluid_select == 23) ||
        (fluid_select == 26) || (fluid_select == 19) ||
        (fluid_select == 27) || (fluid_select == 28) ||
        (fluid_select == 29) || (fluid_select == 40) ||
        (fluid_select == 41) || (fluid_select == 42) ||
        (fluid_select == 43) || (fluid_select == 44) ||
        (fluid_select == 47) || (fluid_select == 147) ||
        (fluid_select == 51) || (fluid_select == 60) ||
        (fluid_select == 61) || (fluid_select == 62) ||
        (fluid_select == 70) || (fluid_select == 90) ||
        (fluid_select == 160)|| 
        (fluid_select == 161) || (fluid_select==170));
}

logical solid_body()
{
return ((fluid_select==70)||
        ((fluid_select>=80)&&(fluid_select<=83))||
        (freetest_flag));
/* 70-periodic test  80-Zalesak's problem */
}

/* below are the routines gdd,gst,guu,gvv,gv,gu,gr,grv,gvt,gut,gru,gd.
   These routines read the data in u,v,d1,d2,d3,ut,vt .... given the
   values ix,iy.   If ix,iy fall outside the boundary an extrapolation step is
   done using the boundary type (mx,my) to compute the value to be returned. */

logical at_vertex(ix)
int ix;
{
  if (
      ((use_symmetric)||(use_axisymmetric)) &&
      (ix<=1)  /* test at 1 not 0 due to corner procedure */
     )
   return TRUE;
  else
   return FALSE;
}   /* at_vertex */

logical use_reflect(ix,xflag)
int ix;
logical xflag;
{
  if (xflag)
   return (at_vertex(ix));
  else
   return FALSE;
}

logical use_extrapolate(ix,xflag)
{
  if (xflag)
   {
    if (at_vertex(ix)) return FALSE;
    else if ((mx_type==REFLECT_BDRY)||(mx_type==SYMMETRIC)||
             (mx_type==FREE_BDRY))
     return TRUE;
    else
     return FALSE;
   }
  else
   {
    if ((my_type==REFLECT_BDRY)||(my_type==SYMMETRIC)||
        (my_type==FREE_BDRY))
     return TRUE;
    else
     return FALSE;
   }
}  /* use_extrapolate */

double get_extrapolate(dd,ix,iy,orientation)
double dd[imax][jmax];
int ix,iy;
int orientation;
{
int ix1,iy1,ix2,iy2,cnt;
double dxval,dyval;

  if (orientation<2)
   {
    if (orientation==0)
     { ix1=1;ix2=2;cnt=1-ix; }
    else
     { ix1=mx;ix2=mx-1;cnt=ix-mx; }

    dxval=dd[ix2][iy]-dd[ix1][iy];
    if (iy<=1)
     dyval=dd[ix1][iy]-dd[ix1][iy+1];
    else if (iy>=my)
     dyval=dd[ix1][iy]-dd[ix1][iy-1];
    else
     dyval=(dd[ix1][iy+1]-dd[ix1][iy-1])/2.0;

    if (fabs(dxval)<fabs(dyval))
     return dd[ix1][iy];
    else
     return dd[ix1][iy]-cnt*dxval;
   }
  else if (orientation>1)
   {
    if (orientation==2)
     { iy1=1;iy2=2;cnt=1-iy; }
    else
     { iy1=my;iy2=my-1;cnt=iy-my; }

    dyval=dd[ix][iy2]-dd[ix][iy1];
    if (ix<=1)
     dxval=dd[ix][iy1]-dd[ix+1][iy1];
    else if (ix>=mx)
     dxval=dd[ix][iy1]-dd[ix-1][iy1];
    else
     dxval=(dd[ix+1][iy1]-dd[ix-1][iy1])/2.0;

    if (fabs(dyval)<fabs(dxval))
     return dd[ix][iy1];
    else
     return dd[ix][iy1]-cnt*dyval;
   }
}   /* get_extrapolate */
  
double gdd(dd,ix,iy,distanceflag)
double dd[imax][jmax];
int ix,iy;
logical distanceflag;
{
double m,b,extval;
int ix1,iy1;

  if (point_ok)
   return dd[ix][iy];
  else if (pointx_ok)
   {
    if ((distanceflag)&&(poolflag()))
     {
      if (iy<=0)
       {
        if (dd[ix][1]>0.0)
         return -dd[ix][1-iy];
        else
         return get_extrapolate(dd,ix,iy,2);
       }
      else
       {
        if (dd[ix][my]<0.0)
         return -dd[ix][2*my-iy+1];
        else
         return get_extrapolate(dd,ix,iy,3);
       }
     }  /* if poolflag */
    else if (use_reflect(ix,FALSE))
     return ( (iy<=0) ? dd[ix][1-iy] : dd[ix][2*my-iy+1]);
    else if (use_extrapolate(ix,FALSE))
     {
      if (iy<=0)
       return get_extrapolate(dd,ix,iy,2);
      else
       return get_extrapolate(dd,ix,iy,3);
     }
    else if (my_type==PERIODIC)
     return ( (iy<=0) ? dd[ix][my+iy] : dd[ix][iy-my] );
    else if (my_type==INFLOW)
     {
      if (iy<=0)
       {
        m=(dd[ix][2]-dd[ix][1])/h;
        b=dd[ix][1]-m*h/2.0;
        return m*(iy*h-h/2.0)+b;
       }
      else
       {
        m=(dd[ix][my]-dd[ix][my-1])/h;
        b=dd[ix][my-1]-m*h/2.0;
        return m*( (iy-my+1)*h+h/2.0)+b;
       }
     }
    else
     return 0.0;
   }  /*  ix in the square */
  else if (pointy_ok)
   {
    if (use_reflect(ix,TRUE)) 
     return ( (ix<=0) ? dd[1-ix][iy] : dd[2*mx-ix+1][iy]);
    else if (use_extrapolate(ix,TRUE))
     {
      if (ix<=0)
       return get_extrapolate(dd,ix,iy,0);
      else
       return get_extrapolate(dd,ix,iy,1);
     }
    else if (mx_type==PERIODIC)
     return ( (ix<=0) ? dd[mx+ix][iy] : dd[ix-mx][iy] );
    else if (mx_type==INFLOW)
     {
      if (ix<=0)
       {
        m=(dd[2][iy]-dd[1][iy])/h;
        b=dd[1][iy]-m*h/2.0;
        return m*(ix*h-h/2.0)+b;
       }
      else
       {
        m=(dd[mx][iy]-dd[mx-1][iy])/h;
        b=dd[mx-1][iy]-m*h/2.0;
        return m*( (ix-mx+1)*h+h/2.0)+b;
       }
     }
    else
     return 0.0;
   }   /* iy in the square */
  else
   {
    if (ix<=0) 
     { ix1=2;ix=1; }
    else if (ix>mx) 
     { ix1=mx-1;ix=mx; }

    if (iy<=0) 
     { iy1=2;iy=1; }
    else if (iy>my) 
     { iy1=my-1;iy=my; }

    if (mx_type==PERIODIC)
     { ix=( (ix<=1) ? mx : 1 );ix1=ix; } /* remember ix changed from 0 to 1 */
    if (my_type==PERIODIC)
     { iy=( (iy<=1) ? my : 1 );iy1=iy; } /* remember iy changed from 0 to 1 */
 
    extval=2*dd[ix][iy]-dd[ix1][iy1];
 
    if ((distanceflag)&&(poolflag()))
     {
      if (iy<=1)
       {
        if (dd[ix][iy]>0.0)
         return -dd[ix][iy];
        else
         return extval;
       }
      else
       {
        if (dd[ix][iy]<0.0)
         return -dd[ix][iy];
        else
         return extval;
       }
     }   /* poolflag */ 
    else if (at_vertex(ix))
     return dd[ix][iy];
    else if ((mx_type==FREE_BDRY)||(use_extrapolate(ix,FALSE)))
     return extval;
    else
     return dd[ix][iy];
   }
}    /* gdd */

double gddd(dd,ix,iy)
double dd[imax][jmax];
int ix,iy;
{
  return gdd(dd,ix,iy,TRUE);
}
 
double gst_general(st1,ix,iy,bottom_stream_value1,top_stream_value1)
int ix,iy;
double st1[imax][jmax],bottom_stream_value1,top_stream_value1;
{
double m,b;

  if (point1_ok)
   return st1[ix][iy];
  else if (point1x_ok)
   {
    if ((my_type==REFLECT_BDRY) || (my_type==SYMMETRIC))
     return 0.0;
    else if (my_type==FREE_BDRY)
     {
      if (iy<iymin)
       return (st1[ix][iymin]*fary+h*bottom_stream_value1)/(fary+h);
      else
       return (st1[ix][iymax]*fary+h*top_stream_value1)/(fary+h);
     }
    else if (my_type==PERIODIC)
     return ( (iy<iymin) ? st1[ix][iymax] : st1[ix][iymin]);
    else if (my_type==INFLOW)
     {
      if (iy<=0)   /* point of outflow */
       {
        return 0.0;
       }
      else
       {       /* point of inflow iy>=my */
        return 0.0;
       }
     }
    else
     return 0.0;
   }  /*  ix in the square */
  else if (point1y_ok)
   {
    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (at_vertex(ix)) )
     return 0.0;
    else if (mx_type==FREE_BDRY)
     {
      if (ix<ixmin)
       return st1[ixmin][iy]*(farx)/(farx+h);
      else
       return st1[ixmax][iy]*(farx)/(farx+h);
     }
    else if (mx_type==PERIODIC)
     return ( (ix<ixmin) ? st1[ixmax][iy] : st1[ixmin][iy]);
    else if (mx_type==INFLOW)
     {
      if (ix<=0)    /* point of inflow */
       {
        return 0.0;
       }
      else
       {      /* point of outflow */
        return 0.0;
       }
     }
    else
     return 0.0;
   }   /* iy in the square */
  else
   {
    if (ix<=ixmin-1) ix=ixmin;
    else if (ix>ixmax) ix=ixmax;
    if (iy<=iymin-1) iy=iymin;
    else if (iy>iymax) iy=iymax;
  
    if (mx_type==PERIODIC)
     ix=( (ix<ixmin+1) ? ixmax : ixmin); /* ix from ixmin-1 to ixmin above */
    if (my_type==PERIODIC)
     iy=( (iy<iymin+1) ? iymax : iymin); /* iy from iymin-1 to iymin above */

    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (at_vertex(ix)) )
     return 0.0;
    else if ((mx_type==FREE_BDRY)&&(my_type==FREE_BDRY))
     {
      return st1[ix][iy]*(farxy)/(farxy+h);
     }
    else if (mx_type==FREE_BDRY)
     return st1[ix][iy]*farx/(farx+h);
    else if (my_type==FREE_BDRY)
     {
      if (iy<=iymin)
       return (st1[ix][iymin]*fary+h*bottom_stream_value1)/(fary+h);
      else if (iy>=iymax)
       return (st1[ix][iymax]*fary+h*top_stream_value1)/(fary+h);
      else
       printf("both iy and ix are out of range check bug\n");
     }
    else if ((mx_type==PERIODIC)&&(my_type==PERIODIC))
     return st1[ix][iy];
    else
     return 0.0;
   }   /* at a corner */
}    /* gst_general */

double gst1(st1,ix,iy)
int ix,iy;
double st1[imax][jmax];
{
return gst_general(st1,ix,iy,bottom_stream_value,top_stream_value);
}

double gst3(st1,ix,iy)   /* used in projection when estimating h+far effect */
int ix,iy;
double st1[imax][jmax];
{
return gst_general(st1,ix,iy,0.0,0.0);
}

double gst2(st2,ix,iy)  /* use in projection invert step */
int ix,iy;
double st2[imax][jmax];
{
  /* we use different bc in projection step than in extrapolation */
  /* REMEMBER: we project to get s_t not s, so the total inflow and
     outflow for s_t is 0.0 while it is top_stream_value-bottom_s..
     for s.  */

  if ( ((mx_type==FREE_BDRY)&&(! point1x_ok)) ||
       ((my_type==FREE_BDRY)&&(! point1y_ok))
     )
   {
    /* the CHANGE in time of inflow is 0; EXCEPT when called from
       calc_stream dealing with u and not ut  */

    if (iy==0)
     return bottom_stream_value_global;
    else if (iy==my)
     return top_stream_value_global;
    else return 0.0;
   }
  else
   return gst1(st2,ix,iy);

}    /* gst2 */


double gst(ix,iy)
int ix,iy;
{
  return gst1(st,ix,iy);
}
 
double guu(uu,ix,iy)
int ix,iy;
double uu[imax][jmax];
{
double m,b,ext,farval;
int ix1,iy1;

/*  return (-pi/314)*(iy-51); */

  if (point_ok)
   return uu[ix][iy];
  else if (solid_body() )
   return gdd(uu,ix,iy,FALSE);
  else if (pointx_ok)
   {
    if ((my_type==REFLECT_BDRY) || (my_type==SYMMETRIC))
     return ( (iy<=0) ? uu[ix][1-iy] : uu[ix][2*my-iy+1]);
    else if (my_type==FREE_BDRY)
     {
      if (iy<=0)
       return uu[ix][1]-(uu[ix][2]-uu[ix][1])*(1-iy);
      else
       return uu[ix][my]-(uu[ix][my-1]-uu[ix][my])*(iy-my);
     }
    else if (my_type==PERIODIC)
     return ( (iy<=0) ? uu[ix][my+iy] : uu[ix][iy-my] );
    else if (my_type==INFLOW)
     {
      if (iy<=0)
       {
        m=(uu[ix][2]-uu[ix][1])/h;
        b=uu[ix][1]-m*h/2.0;
        return m*(iy*h-h/2.0)+b;
       }
      else
       {
        m=(uu[ix][my]-uu[ix][my-1])/h;
        b=uu[ix][my-1]-m*h/2.0;
        return m*( (iy-my+1)*h+h/2.0)+b;
       }
     }
    else
     return 0.0;
   }  /*  ix in the square */
  else if (pointy_ok)
   {
    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (at_vertex(ix)) )
     return ( (ix<=0) ? -uu[1-ix][iy] : -uu[2*mx-ix+1][iy]);
    else if (mx_type==FREE_BDRY)
     {
      if (ix<=0)
       return uu[1][iy]*(farx+(ix-0.5)*h)/(farx+0.5*h);
      else
       return uu[mx][iy]/(farx-(ix-mx-0.5)*h)/(farx+0.5*h);
     }
    else if (mx_type==PERIODIC)
     return ( (ix<=0) ? uu[mx+ix][iy] : uu[ix-mx][iy] );
    else if (mx_type==INFLOW)
     {
      if (ix<=0)
       {
        m=(uu[2][iy]-uu[1][iy])/h;
        b=uu[1][iy]-m*h/2.0;
        return m*(ix*h-h/2.0)+b;
       }
      else
       {
        m=(uu[mx][iy]-uu[mx-1][iy])/h;
        b=uu[mx-1][iy]-m*h/2.0;
        return m*( (ix-mx+1)*h+h/2.0)+b;
       }
     }
    else
     return 0.0;
   }   /* iy in the square */
  else
   {
    if (ix<=0) 
     { ix1=2;ix=1; }
    else if (ix>mx) 
     { ix1=mx-1;ix=mx; }
    if (iy<=0) 
     { iy1=2;iy=1; }
    else if (iy>my) 
     { iy1=my-1;iy=my; }
    
    if (mx_type==PERIODIC)
     { ix=( (ix<=1) ? mx : 1 );ix1=ix; }
    if (my_type==PERIODIC)
     { iy=( (iy<=1) ? my : 1 );iy1=iy; }

    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (my_type==REFLECT_BDRY) || (my_type==SYMMETRIC) ||
        (at_vertex(ix)) )
     return -uu[ix][iy];
    else if (mx_type==FREE_BDRY)
     {
      if ((mx_type==FREE_BDRY)&&(my_type==FREE_BDRY))
       farval=farxy;
      else if (mx_type==FREE_BDRY)
       farval=farx;
      else
       farval=fary;
      ext=2*uu[ix][iy]-uu[ix][iy1];
      return ext*(farval-0.5*h)/(farval+0.5*h);
     }
    else if (my_type==FREE_BDRY)
     return 2*uu[ix][iy]-uu[ix][iy1];
    else
     return uu[ix][iy];
   }
}    /* guu */

double gvv(vv,ix,iy)
int ix,iy;
double vv[imax][jmax];
{
double m,b,ext,farval;
int ix1,iy1;

/*  return (pi/314)*(ix-51); */

  if (point_ok)
   return vv[ix][iy];
  else if (solid_body() )
   return gdd(vv,ix,iy,FALSE);
  else if (pointx_ok)
   {
    if ((my_type==REFLECT_BDRY) || (my_type==SYMMETRIC))
     return ( (iy<=0) ? -vv[ix][1-iy] : -vv[ix][2*my-iy+1]);
    else if (my_type==FREE_BDRY)
     {
      if (iy<=0)
       return vv[ix][1]*(fary+(iy-0.5)*h)/(fary+0.5*h);
      else
       return vv[ix][my]*(fary-(iy-my-0.5)*h)/(fary+0.5*h);
     }
    else if (my_type==PERIODIC)
     return ( (iy<=0) ? vv[ix][my+iy] : vv[ix][iy-my] );
    else if (my_type==INFLOW)
     {
      if (iy<=0)
       {
        m=(vv[ix][2]-vv[ix][1])/h;
        b=vv[ix][1]-m*h/2.0;
        return m*(iy*h-h/2.0)+b;
       }
      else
       {
        m=(vv[ix][my]-vv[ix][my-1])/h;
        b=vv[ix][my-1]-m*h/2.0;
        return m*( (iy-my+1)*h+h/2.0)+b;
       }
     }
    else
     return 0.0;
   }  /*  ix in the square */
  else if (pointy_ok)
   {
    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (at_vertex(ix)) )
     return ( (ix<=0) ? vv[1-ix][iy] : vv[2*mx-ix+1][iy]);
    else if (mx_type==FREE_BDRY)
     {
      if (ix<=0)
       return vv[1][iy]-(vv[2][iy]-vv[1][iy])*(1-ix);
      else
       return vv[mx][iy]-(vv[mx-1][iy]-vv[mx][iy])*(ix-mx);
     }
    else if (mx_type==PERIODIC)
     return ( (ix<=0) ? vv[mx+ix][iy] : vv[ix-mx][iy] );
    else if (mx_type==INFLOW)
     {
      if (ix<=0)
       {
        m=(vv[2][iy]-vv[1][iy])/h;
        b=vv[1][iy]-m*h/2.0;
        return m*(ix*h-h/2.0)+b;
       }
      else
       {
        m=(vv[mx][iy]-vv[mx-1][iy])/h;
        b=vv[mx-1][iy]-m*h/2.0;
        return m*( (ix-mx+1)*h+h/2.0)+b;
       }
     }
    else
     return 0.0;
   }   /* iy in the square */
  else
   {
    if (ix<=0) 
     { ix1=2;ix=1; }
    else if (ix>mx) 
     { ix1=mx-1;ix=mx; }
    if (iy<=0) 
     { iy1=2;iy=1; }
    else if (iy>my) 
     { iy1=my-1;iy=my; }
    
    if (mx_type==PERIODIC)
     { ix=( (ix<=1) ? mx : 1 );ix1=ix; }
    if (my_type==PERIODIC)
     { iy=( (iy<=1) ? my : 1 );iy1=iy; }

    if ((mx_type==REFLECT_BDRY) || (mx_type==SYMMETRIC) ||
        (my_type==REFLECT_BDRY) || (my_type==SYMMETRIC) ||
        (at_vertex(ix)) )
     return -vv[ix][iy];
    else if (my_type==FREE_BDRY)
     {
      if ((mx_type==FREE_BDRY)&&(my_type==FREE_BDRY))
       farval=farxy;
      else if (mx_type==FREE_BDRY)
       farval=farx;
      else
       farval=fary;
      ext=2*vv[ix][iy]-vv[ix1][iy]; 
      return ext*(farval-0.5*h)/(farval+0.5*h);
     }
    else if (mx_type==FREE_BDRY)
     return 2*vv[ix][iy]-vv[ix1][iy];
    else
     return vv[ix][iy];
   }
}    /* gvv */

double gv(ix,iy)
int ix,iy;
{
  return gvv(v,ix,iy);
}


double gr(ix)
int ix;
{
  if (use_axisymmetric)
   {
    if (pointx_ok)
     return x[ix];
    else if (ix<=0)
     return -x[1-ix];
    else 
     return x[mx]+(ix-mx)*h;
   }
  else
   return 1.0;
}

double grv(ix,iy)
int ix,iy;
{
  return gr(ix)*gv(ix,iy);
}

double gvt(ix,iy)
int ix,iy;
{
  return gvv(vt,ix,iy);
}

double gu(ix,iy)
int ix,iy;
{
  return guu(u,ix,iy);
}

double gru(ix,iy)
int ix,iy;
{
  return gr(ix)*gu(ix,iy);
}

double gut(ix,iy)
int ix,iy;
{
  return guu(ut,ix,iy);
}

double gd(dnum,ix,iy)
int dnum,ix,iy;
{
  if (dnum==11)    /* unfiltered version used in conserve.c */
   return gddd(d1,ix,iy);
  else if (dnum==1)
   {
    if (brute_force_flag)
     return gddd(d2,ix,iy);   /* current "distance" version of d1 */
    else
     return gddd(d1,ix,iy);
   }
  else if (dnum==2)
   return gddd(d2,ix,iy);
  else if (dnum==3)
   return gddd(d3,ix,iy);
  else if (dnum==4)
   return gddd(d4,ix,iy);
  else return 0.0;
}

double gdt(dnum,ix,iy)
int dnum,ix,iy;
{
double dhold[3],temp_cutoff;
int i;

  if ((use_bmo)||(! three_phase()) )
   return gd(dnum,ix,iy);
  else
   {
    temp_cutoff=density_spread*h;  /* may be different for viscosity */

    for (i=0;i<3;i++)
     dhold[i]=gd(i+1,ix,iy);
    handle_triple(&(dhold[0]),&(dhold[1]),&(dhold[2]),temp_cutoff);
    return dhold[dnum-1];
   }
}   /* function gdt */
    
